Transcriptomic and metabolomic data of goat ovarian and uterine tissues during sexual maturation

The ovaries and uterus are crucial reproductive organs in mammals, and their coordinated development ensures the normal development of sexual maturity and reproductive capacity. This study aimed to comprehensively capture the different physiological stages of the goat’s sexual maturation by selecting four specific time points. We collected samples of ovarian and uterine tissues from five female Jining Gray goats at each time point: after birth (D1), 2-month-old (M2), 4-month-old (M4), and 6-month-old (M6). By combining transcriptomic sequencing of 40 samples (including rRNA-depleted RNA-seq libraries with 3607.8 million reads and miRNA-seq libraries with 444.0 million reads) and metabolomics analysis, we investigated the transcriptomic mechanisms involved in reproductive regulation in the ovary and uterus during sexual maturation, as well as the changes in metabolites and their functional potential. Additionally, we analyzed blood hormone indices and uterine tissue sections to examine temporal changes. These datasets will provide a valuable reference for the reproductive regulation of the ovary and uterus, as well as the regulation of metabolites during sexual maturation in goats.


Background & Summary
Jining Gray (JG) Goat is a local breed found in the southwestern region of Shandong Province, China.It is known for its early sexual maturity, year-round estrus, and high reproductive capacity 1 .In comparison to other breeds, JG goats reach puberty as early as 2 months of age, with sexual maturity occurring significantly earlier (around 3-4 months) 2 .Therefore, major reproductive features, namely ovarian function and hormonal modulation, are already evident during the early growth stages of JG goats.Consequently, these goats can be considered an optimal model for ovarian development examination in livestock at sexual maturity.Furthermore, their exceptional reproductive traits offer more opportunities for studying sexual maturity.
Sexual maturity refers to the stage after birth when an animal undergoes the growth period and reaches the point where it is capable of normal reproduction.The hypothalamus secretes gonadotropin-releasing hormone (GnRH), which induces the pituitary gland to release follicle-stimulating hormone (FSH) and luteinizing hormone (LH), which subsequently accelerates the onset of development, maturation, and ovarian egg release and hormone secretion such as estrogen and progesterone (PROG) 3 .As estrogen levels increase, the endometrium undergoes proliferation and differentiation, entering the proliferative phase.Subsequently, PROG levels rise, leading to the endometrium entering the secretory phase 4 .These changes signify the onset of the menstrual cycle and the arrival of sexual maturity.The ovaries and uterus are vital organs in the goat's reproductive system.Ovaries synthesize and secrete estrogen and progesterone, as well as maintaining fertility through follicular development and ovulation 5 .The uterus is a target organ for the direct action of ovarian hormones and strictly regulates embryo implantation, pregnancy recognition, and the survival and development of embryos.The physiological functions of the uterus greatly impact the reproductive performance of female livestock, including estrus, mating, conception, and embryo development 6 .The synergistic work of the ovary and uterus facilitates a physiological balance and smooth progression of sexual maturation.This, in turn, lays the foundation for the normal development of reproductive capacity in the offspring.Therefore, the synergistic development of the ovaries and uterus is a key process in the sexual maturation and reproductive development of goats.
Genome-wide sequencing is a widely accepted technology used to comprehensively evaluate the simultaneous changes in animals in response to environmental and dietary induced transcriptional alterations.RNA sequencing (RNA-seq) is a powerful technique for identifying differentially expressed genes and novel transcripts in mammalian reproductive tissues.Specifically, this technology has demonstrated its efficacy in pig gonad 7 , bovine granulosa cell 8 ; goat ovary 9 ; sheep ovary 10 and uterus 11 .On the other hand, metabolomics non-specifically identifies and quantitates all low-molecular-weight metabolic end-products (metabolites).It can enhance our comprehension of the downstream metabolic alterations instigated by post-transcriptional regulation, thereby pinpointing the final stage in a series of modifications triggered by external stimuli 12 .Recent studies have shown that oocyte and gonadal development is strictly modulated by an intricate network of metabolic [13][14][15] .Consequently, the integration of whole-genome sequencing technology and metabolomics offers a robust method for studying the intricate and complex interactions between transcriptional regulation and metabolic processes.This integrated approach holds promise for uncovering the molecular mechanisms that influence reproductive development and exploring metabolic pathways associated with biosynthesis.However, there is currently limited research on the coordinated development of goat ovarian and uterine tissues using both whole-transcriptome and metabolomic studies, especially during the specific stage of sexual maturity.Therefore, in this experiment, we collected ovarian and uterine samples from all the animals on specific days: after birth (D1), as well as 2 (M2), 4 (M4), and 6 months (M6) post birth.These sampling points were chosen to represent distinct physiological phases of goat's sexual maturation.Through a comprehensive analysis of hormone indices, tissue sections, metabolomic, and transcriptomic data, our goal was to investigate the reproductive regulation of the ovaries and the uterus, and understand the metabolic regulation mechanisms during the sexual maturation process of JG goats.
We presented data from transcriptome sequencing and metabolomics assessment of JG goats ovarian and uterine tissues during their sexual maturation phase.A total of 40 samples were used, resulting in sequencing data of 3607.8 million reads for rRNA-depleted RNA-seq libraries and 444.0 million reads for miRNA-seq libraries.The data provided allow an evaluation of varying developmental stages, from birth to post-sexual maturity, to explore the changes in gene transcription activity and metabolism over time in the ovarian and uterine tissues of goats.Both raw and processed data are freely available and potentially contributing to the understanding of the dynamic molecular regulation processes during sexual maturation in JG goats.experimental animal and tissue sample collection.The experiment was conducted at the Shandong Jiaxiang JG Goat Breeding Farm (Jining, China).A total of 20 healthy female JG goats were chosen for this study (Table 1), and they were separated into the following 4 age cohorts: D1 (2.60 ± 1.52 days), M2 (2.07 ± 0.04 months), M4 (4.05 ± 0.05 months), and M6 (6.06 ± 0.06 months).Individual age cohort had total of five goats.Food was freely available to all goats, with breeding and management followed the same protocol.Each goat was uniformly slaughtered on its specific cut-off date.Under sterile conditions, the ovarian (O) and uterine (U) tissues from each individual were collected via surgical instruments and rinsed with cold phosphate-buffered saline (PBS).Forty tissue samples were promptly frozen in liquid nitrogen prior to storage at −80 °C until RNA extraction for NGS library construction and metabolite extraction.Additionally, we collected uterine tissues from each age group, preserving them in a 4% formaldehyde fixation solution with the temperature consistently maintained at 4 °C until ready for histological analysis.Where possible, an initial estimate of developmental stage was obtained through dissection and macroscopic examination of the uterus.

Blood sample acquisition and sex hormone content determination. Jugular vein blood (10 mL)
was extracted from the experimental goat and transferred to a non-anticoagulant tube, which was maintained in a 37 °C water bath for 1 hour, prior to a 10-min centrifugation at 3,000 r/min.The resulting supernatant was allocated into 2 mL RNase-free tubes, and instantly frozen in liquid nitrogen, prior to transport to the laboratory for storage at −80 °C with proper labeling (ID and sample category) for measurement of sex hormone concentrations.To ensure accurate quantification of hormone concentrations, we selected high-quality enzyme-linked immunosorbent assay (ELISA) kits from Qingdao Mdbio Biotech Co., Ltd., headquartered in Qingdao, China.These kits were specifically chosen for the accurate measurement of specific goat hormones, namely, GnRH, FSH, LH, estradiol (E2), PROG, oxytocin (OT), prolactin (PRL), and relaxin (RLN), following the specific protocols provided with each kit.All kits employ the sandwich ELISA method to quantify hormone levels 16,17 , calculating hormone concentrations from the optical density (OD values) obtained at 450 nm wavelength through a standard curve.Each sample was measured thrice, and comparison assessment was conducted using the mean value along with its standard deviation.Samples underwent a 5-fold dilution to ensure the measurements fell within the linear range of the standard curve, with all samples exhibiting a linear regression correlation coefficient (R-value) exceeding 0.95, thus ensuring precision and reliability of our obtained data.The results of each hormone detection and the standard curves are detailed in Table 2 and Supplementary Fig. 1.

Histological assessment.
Fixated uterus tissue underwent dehydration and fixation via rising EtOH (70-100%) and xylene concentrations.individual steps were maintained for a minimum of 30 min.After dehydration, tissues were sequentially placed into three paraffin wax baths at 65 °C for at least 1 hour per bath.Subsequently, the wax-soaked tissue is embedded in the JB-P5 embedding machine.The uterus tissues were sectioned using a RM2016 pathology slicer, with a 4μm slice thickness.Dried sections received Hematoxylin-Eosin staining and observation utilized a Nikon Eclipse E100 light microscope, and images capture employed a NIKON DS-U3 Imaging system.We selected 3 discontinuous sections per sample for observation under a light microscope at 20X, 200X and 400X.The sectional results of uterine tissue from four stages of development, as depicted in Fig. 1., unveil critical indicators of uterine growth: endometrial thickness, glands development, and myometrial thickness 18 .From birth to sexual maturity, we observed notable thickening of the endometrium, enhanced structural integrity, increased gland count, and thickening of the myometrium.These transformations demonstrate the uterus's continuous maturation and refinement, aiding in our further evaluation of sexual development in JG goats during sexual maturation.

RNA isolation and quality analysis.
Total RNA isolation from 20 uterine tissue and 20 ovarian tissue samples employed TRIzol (Invitrogen, Carlsbad, CA, USA) 19,20 following the associated directions.TRIzol reagent contains phenol and additives such as 8-hydroxyquinoline, guanidinium thiocyanate, and β-mercaptoethanol, effectively lysing cells and tissues, releasing various types of RNA including non-coding RNA 21 , and inhibiting RNase activity 22 .The experiment was rigorously conducted according to this procedure: (1) For the Trizol-based method, tissue samples of liquid nitrogen homogenate (50-100 mg) were introduced to 1.5 ml Trizol (Invitrogen, California, USA), prior to gentle 5-8 inversions to mix the suspension, and a 5-min maintenance at room temperature (RT) to ensure complete lysis.(2) Add 300 µl chloroform (J.T. Baker, Pennsylvania, USA), then again inverted for 15 secs, prior to a 5-min incubation at RT.Following a 10-min centrifugation at 12,000 × g at 4 °C, the supernatant will divide into three layers: a bottom red phenol-chloroform organic phase, an interphase, and a top colorless aqueous phase, with RNA being predominantly in the aqueous phase.(3) Transfer the aqueous phase to a fresh tube, prior to introducing equal chloroform volume, and repetition of step 2. (4) The top aqueous phase was transferred to a fresh tube containing 500 µl isopropanol (J.T. Baker, Pennsylvania, USA), followed by a 10-min incubation at RT, prior to a 15-min centrifugation at 4 °C, 12,000 × g.Subsequently, the supernatant is discarded.( 5) Precipitate is rinsed in 1 ml 75% ethanol (J.T. Baker, Pennsylvania, USA).Centrifuge at 4 °C, 7500 × g for 5 min, and supernatant is removed.( 6) Precipitate is air-dried at RT for 5-10 mins, then dissolve the RNA precipitate using an appropriate amount of DEPC-treated water (Qiagen, Hilden, Germany).The RNA concentration, purity and integrity were calculated and checked by a NanoDrop 2000 spectrophotometer (Thermo Scientific, Wilmington, DE, USA) and an Agilent 5400 Bioanalyzer (Agilent, Santa Clara, CA).Currently, the RNA Integrity Number (RIN) is widely used to assess the quality of RNA 23 .RIN values classify RNA samples into 10 predefined integrity categories based on the calculated RIN number for each RNA profile, ranging between 1-10.A RIN = 1 represents completely degraded RNA samples and RIN = 10 represents intact RNA sample.Higher RIN values generally indicate better RNA integrity.RIN values ≥ 8 are optimal for RNA analysis (Table 3) 24 .All RNA sequencing was performed using a single sample.
rRNA-depleted RNA-seq and miRNA-seq libraries construction and sequencing.An Epicenter Ribo-Zero ™ Removal Kit (Epicenter, Madison, WI, USA) was employed for rRNA elimination, and subsequent rRNA-free residues were purified via ethanol precipitation.Sequences that met our strict quality standards were utilized in library generation and sequencing.The lncRNA and mRNA libraries were created with 3 μg total RNA and a NEBNext ®  which yielded the final library with strand specificity.After completing library generation, quantification was completed via a Qubit 2.0 Fluorometer (Life Technologies, CA, USA) and diluted libraries to a concentration of 1 ng/ul.Secondly, an Agilent 2100 bioanalyzer (Agilent Technologies, USA) detected the insert size of the library, which was found to be distributed approximately between 250-300 bp.Finally, the qPCR method was used for the precise determination of optimal library concentrations using a Quantification Kit-Illumina NGS Universal (KAPA, # KK4824) on CFX96 Touch Real-Time PCR Detection System (Bio-Rad laboratories, Hercules, CA, USA), ensuring that the effective library concentration was greater than 2 nM.Ultimately, the suitable libraries underwent sequencing on Illumina NovaSeq 6000 platform (Illumina, San Diego, CA, USA) using the PE150 (paired-end 150 bp) strategy.
Forty small RNA libraries were generated via the NEB Next ® Multiplex Small RNA Library Prep Set for Illumina ® (NEB E7300L) as per the associated directions.First, the 3′ and 5′ adaptors were ligated to the 2 µg of total RNA by T4 RNA ligase for each sample.Subsequently, the first strand cDNA synthesis was performed using M-MuLV Reverse Transcriptase (RNase H-) with the adaptor-ligated RNA as a template.The cDNA was then amplified for 18 cycles using LongAmp Taq 2X Master Mix, SR primer for illumina, and index (X) primer.The resulting PCR products underwent an 8% polyacrylamide gel (100 V, 80 min)-based purification.DNA fragments corresponding to 140~160 bp (the length of small noncoding RNA plus the 3′ and 5′ adaptors) were recovered and dissolved in 8 μL elution buffer for miRNA sequencing library construction.Once the library was constructed, the Qubit2.0Fluorometer (Life Technologies, CA, USA) was used for initial quantification.The library was diluted to 1 ng/µl based on the quantitative results.The insert size of the libraries was then detected using an Agilent 2100 Bioanalyzer (Agilent Technologies, USA).The libraries with insert sizes between 18 and 40 bp were accurately quantified using the CFX96 Touch Real-Time PCR Detection System (Bio-Rad Laboratories, Hercules, CA, USA), and the libraries with an effective concentration above 2 nM were used for further sequencing.The qualified libraries underwent sequencing on the Illumina NovaSeq.6000 platform (Illumina, San Diego, CA, USA), utilizing the SE50 approach (single-end 50 bp, SE50).The aforementioned sequencing was completed by Novogene Co., Ltd.(Beijing, China).

Sequencing data analysis.
The raw imaging data from sequencer was transformed to a sequence file via the CASAVA software (version 1.8.2).This file contained data on both sequence and sequencing quality.To ensure high-quality reads, FastQC software (version 0.11.9) was used to perform quality checks on all samples.For rRNA-depleted RNA-seq libraries, the raw reads filtering utilized the fastp software (version 0.23.1) according to the following parameters: fastp -i in.R1.fq -o out.R1.fq -I in.R2.fq -O out.R2.fq -g -q 5 -u 50 -n 15 -l 150-overlap_diff_ limit 1-overlap_diff_percent_limit 10.Next, clean reads were aligned to the reference genome (Capra hircus, ARS1.2) with Hisat2 (version 2.0.5)(parameters:-phred33-rna-strandness RF-dta-cufflinks-un-conc-gz).Transcript were assembled and quantified were utilized using StringTie software (version 1.3.3b).Due to the use of paired-end (PE) sequencing in rRNA-depleted RNA-seq, transcript profiles levels were normalized to FPKM (Fragments Per Kilobase of transcript per Million mapped reads) to facilitate accurate quantification with RSEM (version 1.3.0).
For small RNA libraries, the raw reads were filtered using Cutadapt software (version 1.16) based on the following criteria: (1) exclude low quality (the bases with a sequencing quality (Q) less than 20 account for more than 30% of the entire read); (2) exclude with 10% or more unknown bases; (3) eliminate unique sequences with a length greater than 30 bp or less than 18 bp; (4) filtering out with harboring ploy-N with 5' adapter contaminants, without 3′ adapter, insert fragments and polyA/T/G/C sequences.Using the Bowtie software (version 1.0.1)filter out repeated sequences and ncRNAs (rRNA, snoRNA, snRNA, tRNA) based on RepeatMasker (https://repeatmasker.org/)and Rfam (https://rfam.org/).Following filtration, the unannotated reads underwent alignment with the goat reference genome by BLAST (version 2.7.1).Finally, the software miREvo (version 1.1) and miRDeep2 (version 2.0.0.7) were integrated to identify miRNAs.Due to the use of single-end (SE) sequencing in miRNA-seq, miRNA expression profiles were normalized to TPM (transcripts per million) with RSEM (version 1.3.0),employing a normalization equation as follows: Normalized expression = (mapped read count / Total reads) * 1,000,000.Principal component analysis (PCA) was conducted on the estimated expression values of the 40 samples using the princomp function (https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/princomp) in the R package.To assess inter-sample batch effect, we conducted relative log expression (RLE) analysis in RUVSeq (version 1.36.0)package.RLE analysis calculates the read count log-ratio to the median count across samples for individual mRNA, miRNA, or lncRNA.

Metabolite extraction and LC-MS analysis.
We conducted non-specific LC-MS profiling sample analysis.Following the previous reported method 25 , the tissues (100 mg) were initially ground in liquid nitrogen before being resuspended in pre-chilled 80% methanol solvent by well vortex.This solvent effectively disrupts cell membranes, promotes cell lysis to release metabolites, and efficiently extracts broad spectrum metabolites, namely, both polar and non-polar substances 26 , while also demonstrating strong protein precipitation capabilities 27 .Comparative studies have shown that methanol outperforms other solvents in capturing diverse metabolites from complex biological matrices 28 .Subsequently, samples underwent a 5-min maintenance on ice, prior to a 20-min centrifugation at 15,000 × g, 4 °C.A portion of supernatant underwent dilution to 53% methanol with LC-MS grade water, then, samples were placed in a fresh tube for a 20-min centrifugation at 15,000 × g, 4 °C.Lastly, supernatant was used to conduct LC-MS/MS system analysis.The apparatus and LC-MS setup for this analysis are described below: a Thermo Scientific Q ExactiveTM HF-X mass spectrometer equipped with a dual-sprayer ESI source was connected to a Vanquish UHPLC system from Thermo Fisher, Germany.This Vanquish UHPLC system comprised essential components, encompassing a vacuum degasser, binary pump, thermostated autosampler, and column oven.To comprehensively cover the metabolome, ionization was conducted in both positive and negative ion modes to maximize the identification of two distinct sets of analytes 29 .The setting of instrument parameters was precisely undertaken: a 3.5 kV spray voltage, 320 °C capillary temperature, 35 psi sheath gas flow rate, and 10 L/min auxiliary gas flow rate.The S-lens RF level was accurately maintained at 60, while the temperature of the auxiliary gas heater was precisely set to 350° for optimal performance.Considering the need for effective retention and separation of medium polarity and non-polar metabolites, the Hypersil Gold column (C18) was utilized for sample analysis.The column set at 40 °C (±1 °C), with 0.2 mL/min flow rate, and 17 min run duration.In positive mode, the mobile phase composition for A was 0.1% formic acid in water, while for B it was methanol.In negative ionization mode, mobile phase A involved a solution of 5 mM ammonium acetate particularly with a pH of 9.0, whereas mobile phase B remained methanol.The gradient elution profile was precisely controlled: an initial 2% B for 1.5 minutes, followed by a linear elevation from 2% to 85% B over 3 minutes.This was succeeded by a gradual rise from 85% to 100%  B over 10 minutes, and thereafter a decline from 100% to 2% B over 10.1 minutes to restore initial conditions.Finally, a 2% B equilibration step was implemented for 12 minutes to stabilize the system.
Data processing and metabolite identification.The UHPLC-MS/MS raw data underwent comprehensive analysis through Compound Discoverer 3.1 developed by Thermo Fisher, which executed peak alignment, peak picking, and quantitation of individual metabolites.Critical parameters encompassed a retention duration tolerance of 0.2 min, actual mass tolerance of 5 ppm, signal intensity tolerance of 30%, signal-to-noise ratio of 3, and a minimum intensity threshold of 100,000.It was attempted to normalize peak intensities to the total spectral intensity, and this normalized value was employed for predicting the molecular formula by amalgamating information from additive ions, molecular ion peaks, and fragment ions.Thereafter, the peaks were meticulously matched against the mzCloud (https://www.mzcloud.org/),mzVault and MassList databases.Compounds

Data Records
The raw rRNA-depleted RNA-seq and small RNA-seq read files for ovarian and uterine tissue have been submitted to the NCBI Sequence Read Archive (SRA) under the project numbers PRJNA1091173 30 and PRJNA1091170 31 , respectively.The raw metabolomics data for ovarian and uterine tissues can be accessed on the MetaboLights database under the accession numbers MTBLS9794 32 and MTBLS9795 33 , respectively.All provided information can be adopted without restrictions.
ethical statement.The research protocol received ethical approval from the review board at Shandong Agricultural University (SDAUA-2023-157).

Fig. 2
Fig. 2 The principal component analysis (PCA) of 40 samples based on estimated expression values.(A) The result of the mRNA sequencing of the ovary; (B) The result of the lncRNA sequencing of the ovary; (C) The result of the miRNA sequencing of the ovary; (D) The result of the mRNA sequencing of the uterus; (E) The result of the lncRNA sequencing of the uterus; (F) The result of the miRNA sequencing of the uterus.

Fig. 3
Fig. 3 The normalized analysis results.Samples in different time periods were displayed in different colors.(A) The violin plots of un-normalized sample RLE for ovarian tissue; (B) Violin plots of normalized sample RLE for ovarian tissue; (C) Violin plots of un-normalized sample RLE for uterine tissue; (D) Violin plots of normalized sample RLE for uterine tissue.

Fig. 4
Fig. 4 The principal component analysis (PCA) of the ovarian and uterine tissue samples datasets acquired in positive and negative ion mode.(A) ovary samples in positive mode; (B) ovary samples in negative mode; (C) uterus samples in positive mode; (D) uterus samples in negative mode.

Table 1 .
Sample information in this study.

Table 2 .
Serum hormones results at different developmental stages.

Table 3 .
The RNA integrity numbers of all samples.

Table 4 .
Raw data, clean data, quality and GC content of 40 miRNA-seq libraries.

Table 5 .
Raw data, clean data, quality and GC content of 40 rRNA-depleted RNA-Seq libraries.